import numpy as np
import matplotlib.pyplot as plt#画图的
import math  #计算幂函数的

from decimal import Decimal
from numpy.linalg import *

def Lagrange():

    line_y = np.zeros(len(line_x))#存储y值得数组

    ans = []#存储结果

    for x in line_x:

        result = 0;#result为一个待求得的函数值
        for i in range(n+1):
            numerator = 1  # 分子
            denominator = 1  # 分母
            for k in range(n+1):
                if i == k :
                    pass
                else:
                    numerator = numerator * (x - x_list[k])
                    denominator = denominator * (x_list[i] - x_list[k])
                    #求ln(x)


            result += (numerator/denominator)*y_list[i]#求Ln(x)


        ans.append(result)

    plt.plot(line_x,ans,label='Lagrange')

    error = get_average_error(ans)

    print("Lagrange 平均误差为 ： " + str(error))

def Vandemonde():
    matrix_x = np.zeros((n+1,n+1))# 范德蒙德矩阵
    line_y = np.zeros(len(line_x))
    matrix_y = y_list
    for i in range(n+1):
        for j in range(n+1):
            matrix_x[i][j] = x_list[i] ** j


    # 解范德蒙德多项式
    ans = solve(matrix_x,matrix_y)

    # 画图
    for i in range(len(line_x)):
        for j in range(n+1):
            line_y[i] += ans[j]*(line_x[i]**j)
    plt.plot(line_x,line_y,label='Vandemonde')

    error = get_average_error(line_y)

    print("Vandemonde 平均误差为 ： " + str(error))

def newton():

    '''说明，由于精度问题，当采样点选择过于密集时，会导致小数溢出，产生“-0.0”这种数字'''
    base_mean = get_base_mean()#计算需要的均差
    ans = []
    for x in line_x:
        result = base_mean[0]#初始值为Pn的第一项，存储所有项的和
        tmp = 1#tmp用于计算每一项的具体值
        for i in range(1,n+1):
            tmp = base_mean[i] #tmp等于均差
            for j in range(i):
                tmp = tmp * (x - x_list[j])

            result += tmp

        ans.append(result)

    plt.plot(line_x, ans, label='newton')

def line_inter():
    '''n取100，插值点取10000比较合适'''
    result = []
    for i in range(n):
        left_b = x_list[i]
        right_b = x_list[i+1]
        help1 = right_b-left_b #x2-x1
        help2 = left_b-right_b #x1-x2
        #取得两个节点之间的区间
        for x in line_x:
            if (x >= left_b and x < right_b) or (x == right_b and i == n-1) :
                tmp1 = ((x - right_b)/help2) * y_list[i]#公式右边
                tmp2 = ((x - left_b)/help1) * y_list[i+1]#公式左边
                final_res = tmp2 + tmp1
                result.append(final_res)

    plt.plot(line_x, result, label='line_inter')

    error = get_average_error(result)

    print("line_inter 平均误差为 ： " + str(error))

def Hermite():
    ans = []
    for i in range(n):
        x_0 = float(x_list[i])

        x_1 = float(x_list[i + 1])

        y_0 = float(y_list[i])

        y_1 = float(y_list[i + 1])

        m_0 = float(d_list[i])

        m_1 = float(d_list[i + 1])
        for x in line_x:

            if (x >= x_0 and x < x_1) or (x == x_1 and i == n-1) :

                help1 = 1 + 2*((x - x_0) / (x_1 - x_0))

                help2 = 1 + 2*((x - x_1) / (x_0 - x_1))

                af0 = float(help1) * ((x-x_1)/(x_0-x_1)) * ((x-x_1)/(x_0-x_1))# (x-x1)/(x0-x1) ^2

                af1 = float(help2) * ((x-x_0)/(x_1-x_0)) * ((x-x_0)/(x_1-x_0)) # (x-x0)/(x1-x0) ^2

                bt0 = (x - x_0) * ( (x - x_1)/(x_0 - x_1) ) * ( (x - x_1)/(x_0 - x_1) )

                bt1 = (x - x_1) * ( ( (x - x_0)/(x_1 - x_0) ) ) * ( ( (x - x_0)/(x_1 - x_0) ) )

                point_y = y_0*af0 + y_1*af1 + m_0 *bt0 + m_1 * bt1

                ans.append(point_y)

    plt.plot(line_x, ans, label='Hermite')

    error = get_average_error(ans)

    print("Hermite 平均误差为 ： " + str(error))

def getDQList(DQList:list):
    #获得均差
    DQTabel = []
    getDQ = lambda column,raw:(DQTabel[column -1 ][raw + 1] - DQTabel[column - 1][raw])/(x_list[column] - x_list[0])


    for i in range(0,n):
        DQTabelColumn = []
        for k in range(0,n-i):
            if(i == 0):
                DQTabelColumn.append(y_list[k])
            else:
                DQTabelColumn.append(getDQ(i,k))
        DQTabel.append(DQTabelColumn)
    for i in range(0,n):
        DQList.append(DQTabel[i][0])
    return

def OmigaX(x, index):
    if index == 0:
        return 1
    result = 1
    for i in range(0,index):
        result *= (x - x_list[i])
    return result

def Newton_Interpolation():
    ans = []
    DifferenceQuotientList = []
    getDQList(DifferenceQuotientList)
    for x in line_x:
        result = 0
        for i in range(0,n):
            result += DifferenceQuotientList[i]*OmigaX(x,i)
        ans.append(result)

    plt.plot(line_x, ans, label='Newton')

    error = get_average_error(ans)

    print("Newton 平均误差为 ： " + str(error))


def get_base_mean():
    matrix_base = np.zeros((n + 1, n + 1))  # 范德蒙德矩阵

    for i in range(n+1):
        matrix_base[0][i] = y_list[i]

    for i in range(1,n+1):#i代表第几阶均差
        for j in range(n+1):#j代表xk的k值
            if i > j:
                matrix_base[i][j] = 0
            else:
                val = (matrix_base[i-1][j]-matrix_base[i-1][j-1]) / (x_list[j] - x_list[j-1])
                final = Decimal(val).quantize(Decimal("0.000"))
                matrix_base[i][j] = final
    #计算均差表

    base = []
    for i in range(n+1):
        base.append(matrix_base[i][i])
    #只取均差表对角线上元素

    return base

def get_average_error(result):

    cnt = 0;
    for i in range(num_insert):
        cnt  += (result[i] - y_list_compare[i])

    error_avg = cnt / num_insert

    return error_avg

if __name__ == '__main__':
    global a,b,c,d,e,f,n,m
    a,b = eval(input("请输入插值区间[a,b],中间用逗号隔开:"))
    c,d,e,f = eval(input("请输入标准函数参数F(x)=c*sindx+e*cosfx中的参数，中间用逗号隔开:"))
    n = int(input("请输入采集点个数n:"))
    m = int(input("请输入实验点个数m:"))
    # a = 0
    # b = 1
    # c = 5
    #
    # d = 30
    #
    # e = 15
    #
    # f = 10
    #
    # n = 100
    # m = 1000

    global num_insert

    num_insert = m

    global x_list,y_list,line_x,d_list,y_list_compare,test_x

    #采样n+1个节点
    x_list = np.linspace(a,b,n+1).tolist()

    #定义函数关系获得y值
    y_list = [c*math.sin(d*i)+e*math.cos(f*i) for i in x_list]

    #求导后的函数值
    d_list = [c*d*math.cos(d*i)-e*f*math.sin(f*i) for i in x_list]


    #获得一个区间内的10000个点
    line_x = np.linspace(a, b, num_insert).tolist()

    #测试平均误差的节点
    test_x = np.linspace(a, b, m).tolist()
    y_list_compare = [c * math.sin(d * i) + e * math.cos(f * i) for i in test_x]


    #绘制原函数
    plt.plot(line_x,[c*math.sin(d*i)+e*math.cos(f*i) for i in line_x],label='sourse')

    # 范德蒙德多项式插值
    Vandemonde()
    # 拉格朗日插值
    Lagrange()
    line_inter()
    Hermite()
    Newton_Interpolation()

    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.show()